home *** CD-ROM | disk | FTP | other *** search
/ Nautilus 1992 July / Nautilus-3-8 / Nautilus-3-8.bin / Tools & Utilities / Techy Stuff / Source ƒ / sky source ƒ / MARS.C < prev    next >
Encoding:
C/C++ Source or Header  |  1992-06-16  |  3.8 KB  |  192 lines
  1. #include "sky.h"
  2.  
  3. extern struct marst
  4. {
  5.     float f[2];
  6.     char c[3];
  7. } marst[];
  8.  
  9. mars()
  10. {
  11.     double pturbl, pturbb, pturbr;
  12.     double lograd;
  13.     double dele, enom, vnom, nd, sl;
  14.     double q0, v0, t0, m0, j0 , s0, u0;
  15.     double lsun, elong, ci, dlong;
  16.     double planp[8];
  17.     struct marst *pp = &marst[0];
  18.     double olong;
  19.     double temp;
  20.  
  21. /*
  22.  *    The arguments nnd coefficients are taken from
  23.  *    Simon Newcomb, Tables of the Heliocentric Motion
  24.  *    of Mars
  25.  *    A.P.A.E. VI, part 4 (1895).
  26.  *
  27.  *    Here are the mean orbital elements.
  28.  */
  29.  
  30.     object = "Mars        ";
  31.     ecc = .09331290 + .000092064*capt - 0.077e-6*capt2;
  32.     incl = 1.850333 - 6.75e-4*capt + 12.61e-6*capt2;
  33.     node = 48.786442 + .770992*capt - 1.39e-6*capt2
  34.      - 5.33e-6*capt3;
  35.     argp = 334.218203 + 1.840758*capt + 1.299e-4*capt2
  36.      - 1.19e-6*capt3;
  37.     mrad = 1.52368840;
  38.     anom = 319.529425 + .5240207666*eday + 1.808e-4*capt2
  39.      + 1.19e-6*capt3;
  40.     motion = 0.5240711638;
  41.  
  42.     incl *= radian;
  43.     node *= radian;
  44.     argp *= radian;
  45.     anom = fmod(anom, 360.)*radian;
  46.     motion *= radian;
  47.  
  48. /*
  49.  *    Conventional mean anomalies of perturbing planets.
  50.  */
  51.  
  52.     q0 = 102.28  + 4.092334429*eday;
  53.     v0 = 212.388 + 1.60211831*eday;
  54.     t0 = 358.415  + .98559696*eday;
  55.     m0 = 319.530 + .52402078*eday;
  56.     j0 = 225.209 + .08307904*eday + 0.332*sin((134.4+38.5*capt)*radian);
  57.     s0 = 175.533 + .03344747*eday - 0.808*sin((134.4+38.5*capt)*radian);
  58.     u0 = 74.188 + 0.0117193*eday;
  59.  
  60.     q0 *= radian;
  61.     v0 *= radian;
  62.     t0 *= radian;
  63.     m0 *= radian;
  64.     j0 *= radian;
  65.     s0 *= radian;
  66.     u0 *= radian;
  67.  
  68.     planp[1] = q0;
  69.     planp[2] = v0;
  70.     planp[3] = t0;
  71.     planp[4] = m0;
  72.     planp[5] = j0;
  73.     planp[6] = s0;
  74.     planp[7] = u0;
  75.  
  76. /*
  77.  *    Computation of long period terms affecting the mean anomaly.
  78.  *        4*mars - 7*earth + 3*venus
  79.  *        3*jupiter - 8*mars + 4*earth
  80.  *        2*jupiter - - 6*mars + 3*earth
  81.  *        2*saturn - 2*mars + earth
  82.  *        jupiter - 2*mars + earth
  83.  *        5*saturn - 2*jupiter
  84.  */
  85.  
  86.     anom = anom
  87.      - (37.05 + 13.50*capt)*radsec
  88.      + 0.606*radsec*sin((212.87+119.051*capt)*radian)
  89.      + 52.490*radsec*sin((47.48+19.771*capt)*radian)
  90.      + 0.319*radsec*sin((116.88+773.444*capt)*radian)
  91.      + 0.130*radsec*sin((74.00+163.00*capt)*radian)
  92.      + 0.009*radsec*sin((325.00+753.67*capt)*radian)
  93.      + 0.280*radsec*sin((300.00+40.8*capt)*radian);
  94.  
  95. /*
  96.  *    Computation of elliptic orbit.
  97.  */
  98.  
  99.     enom = anom + ecc*sin(anom);
  100.     do {
  101.         dele = (anom - enom + ecc * sin(enom)) /
  102.             (1. - ecc*cos(enom));
  103.         enom += dele;
  104.     } while(fabs(dele) > 1.e-8);
  105.     vnom = 2.*atan2(sqrt((1.+ecc)/(1.-ecc))*sin(enom/2.),
  106.             cos(enom/2.));
  107.     rad = mrad*(1. - ecc*cos(enom));
  108.  
  109. /*
  110.  *    Perturbations in longitude.
  111.  */
  112.  
  113.     pturbl = 0.043*sin(2.*anom);
  114.     for(;;){
  115.         if(pp->f[0]==0.){
  116.             pp++;
  117.             break;
  118.         }
  119.         pturbl += pp->f[0]*cos(pp->f[1] + pp->c[0]*m0 + pp->c[1]*planp[pp->c[2]]);
  120.         pp++;
  121.     }
  122.  
  123. /*
  124.  *    Perturbations in latitude.
  125.  */
  126.  
  127.     pturbb = 0.;
  128.     for(;;){
  129.         if(pp->f[0]==0.){
  130.             pp++;
  131.             break;
  132.         }
  133.         pturbb += pp->f[0]*cos(pp->f[1] + pp->c[0]*m0 + pp->c[1]*planp[pp->c[2]]);
  134.         pp++;
  135.     }
  136.  
  137. /*
  138.  *    Perturbations in log radius vector.
  139.  */
  140.  
  141.     pturbr = 0.;
  142.     for(;;){
  143.         if(pp->f[0]==0.){
  144.             pp++;
  145.             break;
  146.         }
  147.         pturbr += pp->f[0]*cos(pp->f[1] + pp->c[0]*m0 + pp->c[1]*planp[pp->c[2]]);
  148.         pp++;
  149.     }
  150.     pturbr *= 1.e-6;
  151.  
  152. /*
  153.  *    reduce to the ecliptic
  154.  */
  155.  
  156.     olong = vnom + argp + pturbl*radsec;
  157.     nd = olong - node;
  158.     lambda = node + atan2(sin(nd)*cos(incl), cos(nd));
  159.  
  160.     sl = sin(incl)*sin(nd);
  161.     beta = atan2(sl, sqrt(1.-sl*sl)) + pturbb*radsec;
  162.  
  163.     lograd = pturbr*2.30258509;
  164.     rad *= 1. + lograd;
  165.  
  166. /*
  167.  *    Compute motion for planetary aberration.
  168.  */
  169.  
  170.     temp = motion*mrad*mrad*sqrt(1.-ecc*ecc)/(rad*rad);
  171.     ldot = temp*sin(2.*(lambda-node))/sin(2.*(olong-node));
  172.     bdot = temp*sin(incl)*cos(lambda-node);
  173.     rdot = motion*mrad*ecc*sin(olong-argp)/sqrt(1.-ecc*ecc);
  174.  
  175. /*
  176.  *    Compute magnitude.
  177.  */
  178.  
  179.     lsun = 99.696678 + 0.9856473354*eday;
  180.     lsun *= radian;
  181.     elong = lambda - lsun;
  182.     ci = (rad - cos(elong))/sqrt(1. + rad*rad - 2.*rad*cos(elong));
  183.     dlong = atan2(sqrt(1.-ci*ci), ci)/radian;
  184.     mag = -1.30 + .01486*dlong;
  185.  
  186.     semi = 4.68;
  187.  
  188.     helio();
  189.     geo();
  190.  
  191. }
  192.